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We compute ground and excited state properties of small helium clusters 4 Hejv containing a single 
benzene impurity molecule. Ground-state structures and energies are obtained for TV = 1,2,3,14 
from importance-sampled, rigid-body diffusion Monte Carlo (DMC). Excited state energies due to 
helium vibrational motion near the molecule surface are evaluated using the projection operator, 
imaginary time spectral evolution (POITSE) method. We find excitation energies of up to ~ 23 K 
above the ground state. These states all possess vibrational character of helium atoms in a highly 
anisotropic potential due to the aromatic molecule, and can be categorized in terms of localized 
and collective vibrational modes. These results appear to provide precursors for a transition from 
localized to collective helium excitations at molecular nanosubstrates of increasing size. We discuss 
the implications of these results for analysis of anomalous spectral features in recent spectroscopic 
studies of large aromatic molecules in helium clusters. 
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I. INTRODUCTION 

Helium droplets provide a unique, ultra-cold nanolabo- 
ratory for investigation of a variety of physical and chem- 
ical phenomenal This has been increasingly used in re- 
cent years to analyze the behavior of a wide variety of 
atomic and molecular species in a quantum liquid envi- 
ronment, using spectroscopic techniques to probe both 
the molecular and helium dynamics* » Electronic spec- 
troscopy in particular allows one to probe microscopic 
details of the helium dynamics. In the large droplet 
regime (N > 10 3 ), the laser-induced fluorescence (LIF) 
spectra of molecules in cold helium droplets are usually 
characterized by a sharp zero-phonon line (ZPL) due to 
the transition for the electronic origin, accompanied by a 
broad, phonon wing sideband. The zero-phonon line can 
contain fine structure due to internal molecular transi- 
tions, while the phonon wing structure reflects the cou- 
pling to collective excitations of the surrounding helium. 
Thus, electronic excitation spectrum of glyoxal, a small 6- 
atom molecule (C2H2O2), in 4 He droplets at T — 0.38 K, 
exhibits a distinct phonon wing feature that is separated 
by a gap from the zero-phonon transition. 3 In contrast, in 
3 He droplets the corresponding spectrum shows a phonon 
wing feature but no gap^ The phonon wing structure for 
glyoxal in has been successfully interpreted in terms of ex- 
citation of the collective phonon-roton modes in large 4 He 
droplets p while the lack of a gap between zero phonon 
and phonon wing features in 3 He droplets has been in- 
terpreted as consistent with the presence of particle-hole 
excitations in normal 3 He4> Recent spectroscopic exper- 
iments involving larger organic molecules in 4 Hej\r clus- 
ters have revealed interesting additional features beyond 
these basic elements £ For a number of the larger organic 
molecules studied so far, both the ZPL and phonon wing 
portions of the LIF spectra exhibit sharp peaks superim- 
posed on the underlying features. These additional peaks 
are not compatible with this basic picture of molecular 



electronic excitation coupling to either internal molecular 
transitions or bulk compressional modes of helium p^LSiS 

In addition to these experiments in large helium 
droplets, a new class of size-selective experiments involv- 
ing small numbers (N < 20) of 4 He atoms attached 
to large planar aromatic molecules has also recently 
emergediiSiii These small cluster studies allow one to 
directly observe the size evolution of excited states in- 
volving helium motion, at sizes less than a full solvation 
shell around the molecule. For these small clusters, which 
are better described as weakly bound complexes with he- 
lium than as quantum solvated molecules, a very different 
situation pertains. Starting with N = 1, the experimen- 
tal spectra show discrete lines that generally increase in 
complexity with increasing N, with the number of ob- 
served lines reaching a maximum at N ~ 4 — 5. After 
this, many of the discrete features observed for smaller 
N disappear, until only a few lines persist at N ~ 10. 

The larger organic molecules that have been stud- 
ied experimentally vary in complexity from planar aro- 
matics such as tetracene (a fused conjugated system 
of four six-membered carbon rings connected by com- 
mon bonds)£i£ to more complex heterostructures such 
as phthalocyanine 9 - and indole derivatives; 8 The pres- 
ence of aromatic character due to 7r-electron conjuga- 
tion provides a common feature in these systems. Be- 
cause of their geometry and 7r-electron character, planar 
aromatics such as tetracene are particularly interesting 
for study in helium droplets because they can be con- 
sidered as nanoscale precursors to a bulk graphite sur- 
face and their local quantum solvation structure con- 
comitantly as a nanoscale precursor of the adsorption 
behavior of thin helium films on graphite. A consider- 
able body of literature has been accumulated for helium 
films on graphite, 12 i 13 i 14 i 15 . 16 i 17 i 18 i 19 i 2Q i 21 i 22 so that anal- 
ysis of the solvation structure and helium excitations as 
a function of increasing size of polyaromatic molecule of- 
fers the possibility of developing a microscopic under- 
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standing of the evolution from quantum solvation of an 
isolated molecule, to adsorption and film formation at a 
bulk surface in superfluid helium. 

The unusual features recorded in experimental spectra 
for these organic molecules consist of unidentified peaks 
in the phonon wings of vibronic bands^ and anomalous 
splittings of the zero-phonon linesi^i^ It has been spec- 
ulated that both of these types of features may be due to 
some type of excitation of helium atoms that are localized 
at the molecular surfacejZ'S'S but the true origin of both 
of these types of features is unclear. Supporting evidence 
for excitations of localized helium atoms derives from the 
similarity of the energies for phonon wing peaks with 
low- frequency modes of thin films on graphite (4—11 K) 
observed via neutron scattering^Si as well as from the- 
oretical predictions of spatially localized helium atoms 
at an aromatic ring in the benzene molecule^ Path in- 
tegral calculations reveal this to be a true localization 
of helium atoms that are effectively completely removed 
from participation in the superfluid, thereby constituting 
a single "dead" atom in the surrounding solvation layer. 
The characteristics of excitations of such localized he- 
lium atoms around aromatic substrates are expected to 
be very different from the collective helium modes found 
in large dropletsr&SLSSi but to show increasing similarity 
with surface localized modes in helium filmsSiSl as the 
size of the molecular nanosubstrate increases. 

Theoretical analysis of the phenomenology of these he- 
lium excitations near a molecular nanosubstrate is ren- 
dered difficult due both to the lack of accurate interac- 
tion potentials for these larger molecules with helium, 
and to the need for accurate calculation of excited states 
in a very inhomogeneous quantum liquid environment. 
Only a few direct calculations to elucidate the nature 
of these helium excitations in the presence of a large 
molecular impurity have been reported so far, and these 
have been restricted to very small numbers of helium 
atoms (N — i^2)£2*£L2L Calculations involving doped 
clusters require accurate helium-impurity potential en- 
ergy surfaces, of which there is a definite lack for the 
larger organic molecules. Thus, most calculations with 
large molecules resort to the use of simple atom-atom 
pair potential models. Given a potential surface, excited 
state calculations employing basis set methods are still 
limited to small N, where N is typically no greater than 
twopiMiyi For larger numbers of helium atoms that are 
however not yet in the regime of bulk films, a practical 
approach is provided by calculations for excited states 
that are based on quantum Monte Carlo methods. 

In this work, we present a series of such quantum 
Monte Carlo studies of the ground and excited states 
of small helium clusters (N = 1,2,3,14) containing 
a benzene impurity. Benzene is a simple planar aro- 
matic molecule, and it can be viewed as a primitive sp 3 - 
hybridized unit of a bulk graphite surface. Thus, it is 
the logical starting point for the kind of systematic study 
of helium adsorption on molecular nanosubstrates men- 
tioned above. We note that this notion has been pursued 



extensively for the heavier noble gases^ but it is only 
recently that experimental data for helium has become 
available. We employ here both the well-known ground- 
state diffusion Monte Carlo method and the projection 
operator imaginary time spectral evolution (POITSE) 
method for excited states^ The latter approach allows 
exact excitation energies to be obtained, provided a high- 
quality trial function approximating the ground state is 
available. We calculate a range of helium excitation en- 
ergies at T — K and analyze these in terms of their 
spatial nature and extent. The previous path integral 
calculations for the He3g-benzene system at a temper- 
ature of T = 0.625 K have shown that the effect of the 
strength and strong 7r-anisotropy of the He-benzene in- 
teraction serves to localize a single helium atom on each 
side of the benzene surfaced That is, in the Feynman 
path integral representation, a single helium atom at- 
tached to the benzene surface is completely removed from 
permutation exchanges with the surrounding helium en- 
vironment. In this work, we find a set of collective helium 
vibrations which have energies of up to ~ 23 K above the 
ground state, and can be characterized in terms of their 
transformation properties under the symmetry group of 
the He-benzene interaction potential. We show that even 
for sizes as small as N = 14 helium atoms, one can dis- 
tinguish a subset of higher-energy excitations localized 
near the global minimum of the hclium-bcnzcnc interac- 
tion potential, from a subset of lower-energy collective 
excitations which are delocalized around the periphery 
of the molecule surface. The former appear to be asso- 
ciated with the localized helium density identified from 
path integral calculations in Ref. |2^ 

Sec. [H] begins with a discussion of the model Hamil- 
tonian and potential surface. Technical details of the 
quantum Monte Carlo methodology are presented in 
Sec. 11111 Our results for ground and excited states of 
4 Hejv-benzene (N = 1,2,3, 14) are presented in Sec. W\ 
where we analyze the nature of these molecule-induced 
localized states around the benzene impurity and, and in 
Sec. lVII we discuss the implications for helium excitations 
on larger aromatic molecules. 



II. THE MODEL HAMILTONIAN 

We treat the 4 HeAr-benzene cluster as a collection of 
TV indistinguishable helium atoms, and a rigid benzene 
molecule which is free to translate and rotate in space. 
Thus the benzene intramolecular degrees of freedom are 
held fixed, which implicitly assumes a separation between 
the helium motion and the benzene vibrational modes. 
The positions of the N helium atoms in a space-fixed (SF) 
frame of reference are denoted as {Ri, R2, ■ ■ ■ , Rjv}) and 
the SF position of the benzene center-of-mass is denoted 
as R/. The Hamiltonian H is (3 AT -I- 6)-dimensional: 
there are 3N helium translational degrees of freedom, 
plus an additional six dimensions due to translation and 
rigid-body rotation of the benzene molecule. It takes the 
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form 



JJ — y(imp) _|_ J>(He) _|_ y ^ 



(1) 



where T( im P), f^ c \ and V are the impurity rigid-body 
kinetic energy, helium kinetic energy, and total poten- 
tial energy, respectively. The impurity kinetic energy 
is most easily expressed by introducing a moving body- 
fixed (BF) frame (x, y, z), whose origin relative to the SF 
frame is fixed at the benzene center-of-mass R/. The Eu- 
ler angles (ip, d, x) specify the orientation of the BF axes, 
which are set parallel to the benzene principal axes. The 
benzene kinetic energy is thus 



f (imp) = - DlV j _ ,/„ ( JL 

with the prefactors 
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n 2 
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2L 



(3) 



The first term of T^ lmp ^ is the center-of-mass transla- 
tional kinetic energy for a benzene molecule, where we 
use a mass to/ = 78.114 amu. The remaining terms give 
the benzene rigid-body rotational kinetic energy, where 
the prefactors d a = 0.188 cm' 1 and d c = 0.0938 cm" 1 
are the (oblate) symmetric top rotational constants. The 
Laplacian V 2 is taken with respect to the translations of 
the benzene center-of-mass. The angular second deriva- 
tives are taken with respect to rotations of the benzene 
about the BF axes (x,y,z). Similarly, the next term of 
the Hamiltonian is the helium kinetic energy, 



f (Ho) 



Y 
3=1 



Da = 



2m4 



(4) 



where TO4 is the 4 He mass, and V 2 denotes the Laplacian 
taken with respect to the displacement of helium atom j. 

The final term of Eq. is the model potential energy. 
We use an additive, two-body potential, 



V 



N 



N 
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(5) 
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The helium-helium interaction Vne-Ho is given by the 
semi-empirical HFD-B potential of Aziz et al.^ and de- 
pends only on the distance r.y between helium atoms 
The helium-benzene interaction T4j e -i is most con- 
veniently expressed in terms of helium BF frame coordi- 
nates {i"i, r2, . . . , rjv}. It is an analytical fit to ab initio 
MP2 data of Hobza et aLfiH and consists of a sum of 
atom-atom terms: 



6 6 



(6) 



TABLE I: Helium-benzene potential parameters. The carbon- 
carbon and carbon-hydrogen bond lengths are held fixed at 
rcc and tch, respectively, and all bond angles are 120°. 



ec 

m 
mi 

rcc 
rcH 



15.25 K 
3.59 A 

19.14 K 
2.69 A 
1.406 A 
1.08 A 



The indices a and (3 run over the carbon and hydrogen 
atoms of the benzene molecule, respectively. The helium- 
hydrogen interaction Vh is a standard Lennard- Jones 6- 
12 form which depends on the distance rpj between hy- 
drogen atom (3 and helium atom j, 



Va(rpj) = 4e H 



12 



(7) 



The helium-carbon interaction Vc has a modified angle- 
dependent Lennard- Jones 8-14 form, 



V c (r aj ) = 4e' c 







( a 'c \ 14 









(8) 



where e' c — ec cos 2 9, a' c = ac cos 9, and 6 is the spher- 
ical polar angle with respect to a coordinate frame cen- 
tered on the carbon atom a and parallel to the BF frame. 
The potential fit parameters are listed in Table [I] This 
helium-benzene potential possesses two equivalent global 
minima of —94.97 K at z = ±3.27 A along the benzene 
Cg-axis. There are also 12 equivalent secondary min- 
ima of —62.54 K located between neighboring hydrogen 
atoms, six above and six below the benzene plane. The 
global minima are connected to the secondary minima by 
saddle points of —50.77 K, situated approximately above 
and below the C-H bonds. A cut of the potential at 
z = 3.27 A is shown in the top panel of Fig. ^ where 
the (a;, y) positions of the carbon and hydrogen atoms 
on the z = plane are also marked. The lower panel of 
Fig.^shows a cut of the potential along the y — plane, 
which reveals one of the secondary minima at x = 2.84 A, 
z = 2.70 A, and one of the saddle points at x = 1.97 A, 
z = 3.54 A. 

For N = 1, the model Hamiltonian of Eq. (JIJ be- 
longs to the molecular symmetry group D 6 h (M) , which is 
the subgroup consisting of only the feasible elements of 
the complete nuclear permutation and inversion group 
(CNPI). 36 This group is isomorphic to the D§h point 
group, which is valid in the limit of small amplitude he- 
lium motions. For N > 1, the symmetry group is iden- 
tically D 6 h(M) due to Bose symmetry, since the only 
allowed rovibrational states are those obtained from the 
tensor product of irreducible representations of Dg/j(M), 
with irreducible representations of the symmetric group 
Sn which remain unchanged under permutation of iden- 
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FIG. 1: Top: cut of the 4 He-benzene potential through the 
global minima at z = 3.27 A. The carbon and hydrogen 
(x, y) positions are marked with solid circles and connected 
by dashed lines. Bottom: cut of the same potential at y = 0. 
Three representative sites sy , sy/ , sy// are marked, which cor- 
respond to centers of the helium-benzene localization factors 
centered near the global, secondary, and tertiary potential 
minima, respectively fSecs. HVl and lV^ll . 



tical 4 He nuclei, i.e. the totally symmetric representation 
ofS N . 



III. COMPUTATIONAL METHODOLOGY 

In this section we present a summary of the numeri- 
cal methods employed in our study of the He^-benzene 
system. Variational estimates of ground-state energies 
are obtained from variational Monte Carlo (VMC), and 
exact ground-state properties are computed from diffu- 
sion Monte Carlo (DMC) with importance samplingi2L2£i 
The projection operator, imaginary time spectral evolu- 
tion (POITSE) approach is used to obtain excited-state 
energies^ Both the DMC and POITSE methods can 
provide results that are exact, to within a systematic 
time step error. 



A. Ground states: variational Monte Carlo (VMC) 

The starting point for our study of the 4 He7v-benzene 
system begins with the Rayleigh-Ritz variational theo- 
rem, 



En 



(^ T \H\^ T ) 



> E n 



(9) 



where the trial energy Et with respect to some param- 
eterized trial function ^>t represents an upper bound to 
the exact ground-state energy Eq. In the coordinate rep- 
resentation, this becomes 



Ex — 



J dH^ T {K)E L (K) 
JdTZ^ T {TZ) : 



(10) 



where 1Z is a (3iV + 6)-dimensional coordinate denoting 
the positions and orientations of all bodies governed by 
the Hamiltonian H. The quantity El is a local energy, 
defined as 



(11) 



An optimized variational upper bound is obtained by 
varying the parameters of to minimize En?. 

For a realistic V-particle helium system, Monte Carlo 
methods offer a practical means to compute the multi- 
dimensional integral of Eq. I |10| l. In variational Monte 
Carlo, expectation values of observables are evaluated as 
averages over the normalized distribution 



(12) 



which is numerically represented as a discrete ensemble 
of Monte Carlo walkers {IZk}- Typically, we use an en- 
semble size of n w = 2000 walkers. This distribution can 
be generated from a Metropolis walk, in which a move 
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TZ — > TZ' is proposed with a transition probability con- 
sisting of factors having the form 



G D (Rj ^r;,at) 
1 



4ttDAt 



3' 
3/2 



exp ■ 



4DAt 



(13) 



The vector quantity F is commonly referred to as the 
"quantum force" , and is chosen to be 



Fj(K) = 24v 1 (7e)V i * T (ft). 



(14) 



Gd can be viewed as the Green's function associated 
with a diffusion/drift process in the presence of an ex- 
ternal field F, over a time step of At. All VMC- 
based computations reported here use a time step of 
At = 75 — 100 Hartrcc -1 . For the 4 Hejv-benzene sys- 
tem the full (3N + 6) -dimensional transition probability 
is comprised of N such factors with diffusion constant 
D = Da corresponding to the N helium atoms, a fac- 
tor with D = Dj corresponding to the benzene impurity, 
and three analogous one-dimensional factors correspond- 
ing to individual rotations of the benzene about each of 
its three principal axes. A proposed move TZ — > TZ' is sub- 
sequently accepted with probability p = min(l, q), where 
the acceptance ratio q is 



^ 2 T {TZ')G D {TZ' -^TZ,At) 
^ 2 T {TZ) G d (TZ^TZ',At) 



(15) 



With the choice of transition and acceptance probabil- 
ities given by Eqs. (|T^|l - lfT5|) . the Metropolis walk con- 
verges to the asymptotic distribution ^^(TZ). At this 
point, expectation values for quantities such as the vari- 
ational trial energy Et are sampled from this distribution 



/ 1 n w 

E T = (E L ) = ( — ^£ L (ftfc) 



(16) 



k=l 



walk 



where the notation (. . .) wa ik denotes a statistical average 
over the course of the equilibrated Metropolis walk. This 
is performed by sampling the quantity given in (. . .) wa ik 
every M time steps apart, with M chosen to be longer 
than the autocorrelation length in order to minimize cor- 
relation biases in the error estimates; typically M = 50 
for our VMC computations here. 



method derives from the imaginary time (r 
Schrodinger equation, 



d^(TZ, t) 



= (H-E ref )iH(n,T), 



it/h) 



(17) 



where ^{TZ, r) is the many-body wave function, and E re f 
is an arbitrary constant shift in the energy spectrum. 
Importance sampling is introduced by multiplying both 
sides of Eq. (|17|) by a trial function ^>t(TZ), and rewriting 
in terms of f(7Z,r) = ^t(TZ)^(TZ, t ) to obtain a set of 
equations having the form 



df(K,r) 
dr 



+ {E L (K)-E ref ]f(K,T). (18) 



The stationary solution is now the normalized "mixed" 
distribution 



V T (1Z)MK) 
J dTZ^ T (TZ)M^)'' 



(19) 



where tfio is the exact many-body ground-state wave func- 
tion. In DMC, the approximate short-time Green's func- 
tion G which generates this distribution is 

G(1Z -» 1Z', At) » G D (1Z -> 1Z' . At)G b (1Z -> 1Z' . At). 

(20) 

Here, Gd is the diffusion/drift Green's function of 
Eq. H13[l . and Gb has the form 



G B {1Z 



exp 



R',Ar) = 
E L (1Z) 



E L (TZ') 



-E 



rcj 



(21) 



All DMC-based computations here use a time step of 
At = 25 Hartree -1 . A proposed move TZ — » TZ' is ac- 
cepted or rejected according to Eq. Ijl5|l . to ensure that 
the exact ground state is sampled in the limit of per- 
fect importance sampling, i.e. when ^t(TZ) — (f>o(lZ). 
For a non-exact trial function however, the approximate 
Green's function of Eq. I|20[l results in a systematic time 
step error. 

Operationally, the importance-sampled DMC proce- 
dure is similar to that described previously for VMC, 
except now associated with each walker TZk is a cumula- 
tive weight Wk due to the action of Gb' 



J\G B {TZk^K,nAT). 



(22) 



B. Ground states: diffusion Monte Carlo (DMC) 

With diffusion Monte Carlo (DMC), one can move 
beyond the variational level of theory to obtain exact 
ground-state expectation values for the energy and quan- 
tities which commute with the Hamiltonian H. This 



The efficiency of the DMC method can be significantly 
improved by replicating walkers with large Wk, and de- 
stroying walkers with small Wk- At every 10 — 50 time 
steps, the ensemble is checked for walkers whose weight 
exceeds the empirically set bounds w m i n and w max . A 
walker TZk with weight Wk > w max is replicated into 
rik = int(wfc + it) walkers, where u is a uniform random 



number on [0, 1). These new walkers are then assigned 
the weight Wk/n k . A walker lZ k with weight w k < Wmin 
is destroyed with probability 1 — w k \ otherwise its weight 
is set to unity. This is to ensure that the branching 
scheme, on average, does not artificially alter the ensem- 
ble sum of weights Wtot = J2 w k- The parameters w m in 
and Wmax are chosen empirically to maintain a stable 
DMC walk with respect to Wtot and the ensemble size; 
here we use w m i n = 0.1 — 0.3 and w max = 2.0 — 2.2. Ad- 
ditionally, we also vary the reference energy E re t during 
the course of the walk according to 



E, 



W tot {T) 



Wtot(r + AT) 



(23) 



Here, the imaginary-time dependence of the E re f and 
Wtot is explicitly written. We set E re f initially to be the 
starting ensemble average for the local energy (El). The 
empirical update parameter rj is chosen to be as small as 
possible, typically t]/At — 0.0004 — 0.004, which gives 
fluctuations in W to t of roughly 5%. 



computing instead the reweighted average 

Efeli w k (r)w k (T + M'Ar)* T 1 (^ fe )6* T (^ fc ) 



(0) = 



Y,l=iMr)w k {T + M' At) 



f T (R) 



= {4> Q \0\(f)o), 



(27) 



(28) 



(29) 



where M 1 is taken to be as large as possible, typically 
Af = 1500. When branching processes are incorpo- 
rated in the DMC, one needs to take care to keep track 
of walkers that descended from a walker TZ k at time r. 
All helium densities reported in this work (N < 14) are 
computed by reweighting walkers TZ k by their descendant 
weights, i.e. using Eq. i|27|) . However, the statistical noise 
seen in the densities computed in this manner is much 
greater than from densities obtained using the mixed es- 
timator, and thus this approach becomes problematic for 
4 HeAr-benzene at still larger sizes. 



Once the ensemble converges to the mixed distribution 
^!t{TV)4'o{T^)i mixed expectation values for a Hermitian 
observable O can be obtained as 



Excited states: projection operator imaginary 
time spectral evolution (POITSE) 



(O)mix = l^T\ Wk^iTZ^O^TiTZk) ) (24) 

\ Wtot i, / 
\ ft=l I walk 

/ dn^ T {n)(i)Q{n)^ l {n)d^ T {n) 



J dK^ T {K)4> (n) 

6 |6|*t) 



|*T> 



(25) 
(26) 



In our ground-state DMC calculations, this average is 
computed by sampling O every M — 150 time steps 
apart, which is larger than that used for VMC because 
successive DMC iterations are more strongly correlated 
due to the smaller DMC time steps that we use here. 
In the case where <j>o is an eigenstate of O, i.e. when O 
commutes with the Hamiltonian H , this procedure yields 
exact expectation values. All ground-state energies re- 
ported in this work are derived from the mixed estimator 
for the local energy, and are thus exact. For quantities 
which do not commute with H, such as the local helium 
density operator p(r) — J2j <^( r — r .j')> the mixed estima- 
tor is biased by the trial function ^t- But as pointed 
out by Liu et al.^ the asymptotic weight Wk(r — > oo) of 
a walker IZk is proportional to <f) (lZk) /^riTlk)- Thus, 
in principle this trial function bias can be eliminated by 



The POITSE approach is a DMC-based method which 
can provide exact excited state energies, subject to a sys- 
tematic time step biasi^ It begins with the DMC evalu- 
ation of the imaginary-time correlation function 



«(r) 



(* T |ie-^- E -/) T i+|* T ) 



(30) 

|i t |* T )| 2 e- (iS '^ Bo)T + 0(a; 2 ) (31) 



(^ T | e -(«~s„/)r|^ T ) 



where 



(32) 



Here, {4> n } and {E n } represent a complete set of energy 
eigenstates and eigenvalues of the Hamiltonian H, re- 
spectively, and A' is an operator chosen to connect ^>t 
to the excited state(s) of interest <f> n . should be a 
good approximation to the exact ground state <pQ. To 
second order in x, k(t) is a superposition of exponential 
decays whose decay rates correspond to the energy differ- 
ences E n - E 40Al The VMC procedure of Sec- lrrTAl is 
used to generate an initial ensemble of n w = 2000 walkers 
distributed according to f(lZ,0) oc ^^(TZ). This distri- 
bution is then propagated in imaginary time using the 
DMC procedure outlined in Sec. IIII Bl during which the 
correlation function k(t) is sampled from the DMC walk 
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k(t) 



1 



W tot (r) 



n w (r) 

ti{K^)A{KV)w k ,{T). 



(33) 



The index kl denotes walkers at time r which descended 
from a parent walker k at time r = 0, and here we also 
explicitly emphasize the imaginary-time dependence of 
the DMC quantities Wtot, n w , and Wk. 

An inverse Laplace transform of Eq. I|31|l yields the 
spectral function 

k(E) <xJ2\(MA f \yT)\ 2 5(E - E n + Eo), (34) 



whose peak positions give the excited state energies 
E n — Eq. The additive contributions 0(x 2 ) are neglected 
in Eq. I|34fl . since they do not affect the positions of the 
peaks associated with E n — Eq, and in practice, their 
spectral weights are negligible for a reasonable choice of 
\]/ Tl ^M£ The numerical inverse Laplace transform is per- 
formed using the Maximum Entropy Method (MEM), in 
which the Laplace inversion is formulated as a data anal- 
ysis problem in terms of Bayesian probability theory^ 
We adopt here the vector notation k = {n(Ei)AE} to de- 
note the discrete spectral function specified at intervals 
of width AE, and k = {^(rJAr} to denote the correla- 
tion function sampled at intervals of At. The Bayesian 
approach begins with the posterior distribution function 



p(n\k, a, m, I) oa e 



aS-L 



(35) 



The posterior p represents the probability of obtaining 
the image k, given the Monte Carlo data k, the param- 
eter a, an initial guess for the image m (also referred 
to as the default model), and any other relevant back- 
ground information /. We use a constant value for the 
default model m. The quantity 5* is the Shannon-Jaynes 
information-theoretic entropy, 



m log — 



(36) 



Carlo noise. For a > 0, the search for k requires that 
the entropy S be simultaneously maximized while x 2 is 
minimized. Thus, a can be viewed as a regularization 
parameter, which stabilizes the least-squares fit by con- 
straining the minimization of x 2 ■ 

In the ideal situation where the spectrum consists of a 
single, sharp peak, we take the excitation energy E n — Eq 
to be the first moment of the spectrum, 



E n — E — (E) — 



J2i E *Ki 



(38) 



The corresponding estimate for the variance of the mean 
<j then follows from the usual procedure for the propa- 
gation of uncertainties, 



C 



dk 



(39) 



which requires knowledge of the covariance matrix C for 
the image k. By approximating p as a sharply-peaked 
Gaussian in the neighborhood of k, the covariance of the 
image becomes^ 



-(aVVS - VVL)~ 



(40) 



When the spectrum is composed of multiple, well- 
separated peaks, the mean value of the peak and the 
variance in the mean is obtained in the same manner 
as given by Eqs. i|3S|l - l|59"|l . except that the summation is 
taken only over the spectral feature of interest. Note that 
the approach taken here in evaluating error estimates dif- 
fers somewhat from the Gaussian approach advocated in 
Ref.^3 I n the Gaussian approach the variance in E 71 —Eq 
is associated with the peak width, and thus does not nec- 
essarily scale as 1/Nd, where Nd is the number of DMC 
decays used to compute statistics for the covariance ma- 
trix C of Eq. I|37|) . Here, we obtain E n — Eq from the 
mean peak position, and then evaluate the variance of the 
mean according to Eqs. I|38ll - (|39[) . which by construction 
scales as 1/N d . 



and L = x 2 /2 is chosen to be the usual Gaussian likeli- 
hood. In matrix notation, x 2 takes the form 



X 



(Ck - kf ■ CT 1 • (Ck - k), 



(37) 



where C is the Laplace operator, and C is the covariance 
matrix for the data k. In this work, C is obtained from 
averaging Nd — 256 decays, which is sufficient to give 
well-converged spectra for N < 3. A search is then made 
for the image k which maximizes the posterior p, using 
a search algorithm due to Bryan. 43 In the limit a = 0, 
this reduces to a standard least-squares data-fitting pro- 
cedure involving the minimization of x 2 , which is numer- 
ically unstable when one seeks to infer a free-form solu- 
tion for the image from data with non-negligible Monte 



IV. TRIAL FUNCTIONS AND EXCITATION 
OPERATORS 

While in principle a fully converged DMC result should 
be independent of the trial function \I/t, a good choice 
of which closely approximates the exact ground state 
4>q can significantly improve computational efficiency. On 
the other hand, a poor choice of ^>t can produce mislead- 
ing results. In this study, we use analytical forms moti- 
vated by basic physical considerations. The trial function 
is a product of two- body factors, 



JV 



N 



*T(rc)=rufo)nx(r«). 



(41) 
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Here, x describes helium-helium correlations, and has the 
McMillan formic 



exp 



CHc 



(42) 



The function £ describes helium-benzene correlations, 
and is a product of helium-benzene factors defined in the 
BF frame, 



£(ri) = Il£> 



(43) 



7=0 



The first of these factors £o controls the behavior of the 
trial function at short and long hclium-bcnzcnc separa- 
tions: 



exp 



6 

E cc 
„6 



6 

E CH 

1 ' otj /3=1 0j 



ao{x] + y 2 ) - c z 



(44) 



The terms involving the parameters cc , ch ensure that 
goes to zero as a helium atom approaches a carbon or 
hydrogen atom, where r a j and rpj denote helium-carbon 
and helium-hydrogen distances, respectively. The r~ 6 
and r~ 5 forms are chosen to cancel the leading singular- 
ity in the local energyai due to the helium-benzene pair 
potential of Eq. ©. 

With £ = £o alone, it is known that the resulting trial 
function describes a ground state with "liquid" -like char- 
acteristics. If the true ground state has "solid" -like qual- 
ities, i.e. individual helium atoms are strongly localized 
around the molecular impurity, such a state will not so- 
lidify out of a liquid-like trial function at the variational 
level of theory, within reasonable VMC simulation times. 
One simple way to incorporate solid-like features a pri- 
ori is to introduce a set of sites {s 7 } fixed in the BF 
frame. The remaining factors £ 7 (7 > 0) serve to localize 
helium atoms near these sites. In this work we use the 
exponential of a Gaussian to effect this localization: 48 



6y( r i) = CX P ( c 



yl r j- s 



(45) 



The specification of the sites {s 7 } is given later in 
Sec. IV Al A more sophisticated variational approach 
would involve the extension of this trial function to 
shadow functions, in which the {s 7 } are treated as sub- 
sidiary variables, thus obviating the need for an a pri- 
ori specification of {s 7 }. 49 While a shadow trial function 
is expected to yield a significantly better description of 
the ground state, we do not employ this approach here 
because it introduces additional complexity due to the 
need to sample the subsidiary variables {s 7 } in VMC 
and DMC. 

In the POITSE formulation, an initial ansatz for the 
excited state is generated by the action of an excitation 
operator on ^t, and exact excited-state energies are 



then extracted from the DMC imaginary-time evolution 
of this initial state. Due to the inherent numerical diffi- 
culties in the MEM inversion when the target image k(E) 
contains multiple, closely-spaced peaks of comparable 
spectral weight, we seek operators A 1 which connect ^>t 
to only one or a few energetically well-separated states. 
Thus, one consideration is to choose A^ to transform as 
an irreducible representation of the DQh(M) group, such 
that the spectral weight (4> n \A^\^T) is only non-zero for 
states <j) n which transform identically. For systems with 
a relatively high degree of symmetry, such a choice of 
A^ would significantly reduce the number of individual 
decays contributing to k(t). 

In the BF representation, a convenient primitive for 
A' are provided by the regular spherical harmonics: 



-Rzm(r) 



47T 

21 + 1 



1/2 



Yi m {B,<t>). 



(46) 



Similar operators have been used previously to study ro- 
tationally excited states of small 4 Hejv and (H^n clus- 
ters (N = 7, 40) i 50 i 51 For a pure cluster, this yields trial 
excited states which are simultaneous eigenstates of the 
square of the total helium angular momentum L 2 , and 
its SF z-component L z . Addition of a benzene impurity 
lowers the symmetry to DQh(M)- In this case the set of 
regular spherical harmonics can be symmetrized by tak- 
ing linear combinations of Ri m and Ri t - m to obtain real 
operators which transform as irreducible representations 
of D 6h (M): 



iRlr 



l) m -Rim + Rl,-r 



l) m Rlm - Rl,- 



(47) 



(48) 



The subscripts c and s indicate proportionality to 
cos(m0) and sin(m</>), respectively. The various Ri m are 
listed in Table II in terms of BF Cartesian coordinates 
for up to I — - 4, along with their corresponding symme- 
try labels under D 6h (Al). Note that for up to I = 4 
under Deh(M), Ri mc and i?/ ms are doubly degenerate 
except for m — 3, and thus from this point on we drop 
the subscripts c and s for the degenerate functions ex- 
cept when the distinction is necessary. For N > 1, these 
one-particle operators are Bose-symmetrized over the N 
identical bosons to give single-excitation operators: 



JV 



(49) 



Application of AJ on &t then gives the trial excited 



state: 



(50) 

JV N 

= X)£(ri)£(r 2 ) . . . Ri m (rMrj) ■ ■ ■ t(r N ) J] x(r y ). 

(51) 
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TABLE II: Symmetrized regular spherical harmonics Ri m of 
Eq. (1461 - 14811 . in terms of BF Cartesian coordinates (to within 
a constant real factor). The last column gives its correspond- 
ing symmetry V according to the irreducible representations 
of the D$h(M) group. 



Rim 








r 


Rio 




z 




A 2u 


Rue Rus 




x, y 




Eiu 


R20 




x 2 + y 2 - 2z 2 




Aig 


R21C, R21S 




xz, yz 




Elg 


R22C, R22S 




2 2 
x - y , xy 




E 2 g 


Rao 




/ 2 , 2 2 2\ 

z{x +y --sZ) 
-4z 2 ), y(x^ + y 2 




A 2u 


R31C, R31S 


x(x' 2 + y 2 


-4z 2 ) 


Eiu 


R'i2c, R'A2s 




z(x 2 — y 2 ), xyz 




E 2u 


R 3 3c 




x 3 - 3xy 2 




B-2u 


R33s 




3x 2 y - y 3 




Biu 


R40 


8z 2 [z 2 - 


3(x 2 + y 2 )} + 3(x 2 


+ y 2 ) 2 




Raic Ras 


xz[4z 2 — 3(ar 
x 2 -y 2 )(6z 2 


2 + y 2 )], yz[4z 2 -: 


^{x 2 +y 2 )] 
'■~x 2 ~y 2 ) 


Elg 


Rl2c, Ri2s ( 


-x 2 -y 2 ), xy(6z' 


E 2 g 


Ri3c 




z(x 3 - 3xy 2 ) 




Big 


Ri3a 




z(3x 2 y-y 3 ) 




B 2 g 


RlAc, RiAs 


(x 2 -y 2 


) 2 - (2xy) 2 , xy(x 2 


-y 2 ) 


E 2 g 



TABLE III: Parameters for * T of Eqs. l^-lgKJ, in atomic 
units. The position of a site s 7 is listed in terms of its (x, y, z) 
Cartesian coordinates in the BF frame. Only one center from 
a set of sites which are equivalent by symmetry is given, and 
the remaining can be generated by application of the elements 
of the Deh{M) group. 





I 


2 


r> 



14 


CHe 




AACiA 1 


AAf\A 1 




CQ 


6000.0 


6000.0 


6590.9 


8217.7 


CH 


8000.0 


8000.0 


3519.7 


2546.2 


do 


0.05 


0.05 


0.0078624 


0.014378 


CO 


0.06 


0.06 


0.01441 


0.0073448 








2.79 


1.9317 


<Zy 






0.090199 


0.15613 








(0.0,0.0,7.0) 


(0.0,0.0,7.0) 


Cy/ 






1.7103 


1.5165 


(Jy/ 






0.090199 


0.095327 








(6.0,0.0,5.5) 


(6.0,0.0,5.5) 










2.217 


ay// 








0.030121 










(9.5,0.0,0.0) 



Thus, the action of Aj is to impose a trial nodal struc- 
ture on the helium-benzene factor £, generating a trial ex- 
cited state corresponding to a Bose-symmetrized single- 
excitation excited state. As will be discussed further in 
Sec. IV Bl these operators can be further refined to take 
into account the localization of the helium density near 
the molecule surface. In principle this approach can also 
be readily extended to multi-excitation operators. 

V. RESULTS AND DISCUSSION 
A. Ground states 

Ground-state energies and structures for the HeAr- 
benzene cluster (N = 1,2,3,14) were computed using 
the methodologies described in Secs. lIII Al and HTlBI Our 
general strategy for ground states is to first begin with 
a preliminary unbiased DMC calculation with = 1, 
which gives an asymptotic DMC distribution propor- 
tional to the exact Bose ground state (fo. However, such 
an unbiased DMC approach for helium clusters has a 
number of convergence issues associated with it. In the 
absence of a guiding function, the DMC walk will spend 
too much time sampling unimportant regions of configu- 
ration space, and thus unbiased DMC can give an energy 
which is not converged and is higher than the true value 
of Eo. 38 Nevertheless, it can be useful in situations where 
no other information is available, in particular to act as 
a guide for constructing a suitable starting trial func- 
tion ^t- We first compute a reduced ground-state wave 
function 4>o(r) in the BF frame by binning the unbiased 
DMC distribution <fio(lZ) onto a three-dimensional grid. 
Trial function parameters for the helium-benzene factor 



TABLE IV: Ground-state energy per particle Eo/N, in 
Kelvins. The second column from left gives the variational 
energy from VMC with respect to the trial function ^t, pa- 
rameterized by values given in Table |ffl| Unbiased DMC 
refers to DMC with * T = 1. IS-DMC is the best estimate for 
the exact ground-state energy, using the trial function whose 
corresponding VMC energy is given in the second column. 
The value in parenthesis represents the standard deviation in 
the last significant figure. 



N 


VMC 


unbiased DMC 


IS-DMC 


1 


-1.9(2) 


-26.72(7) 


-26.83(4) 


2 


-4.0(2) 


-27.0(1) 


-27.12(3) 


3 


-9.6(2) 


-23.8(1) 


-23.85(2) 


14 


-9.6(3) 


-16.11(5) 


-17.51(1) 



£, Eqs. P%|) - (|4l)|l . are then obtained by fitting to this 
</>o( r )- In some cases, these parameters are further var- 
ied to lower the variational energy, and the resulting trial 
function is then used as input for an importance-sampled 
DMC calculation. The final optimized trial function pa- 
rameters for the various sizes are listed in Table ITTT1 

The VMC and DMC results for the ground-state en- 
ergy per particle Eo/N are summarized in Table ITVl For 
N = 1 — 2, we find that a simple trial function involving 
£ = £0 only, Eq. Ij44(l . is sufficient to give good DMC re- 
sults. With the two global minima occupied at N = 2, 
additional helium atoms begin to occupy the twelve sec- 
ondary potential minima (Fig. ^| . To describe this situ- 
ation for TV = 3, we incorporate localization factors £ 7 , 
Eq. I|45|). centered at sites {sy} and {sy}. The set of 
sites {sy} are centered near the two equivalent global 
minima of the helium-benzene potential, the set {sy} 
are centered near the twelve equivalent secondary poten- 
tial minima. The N = 14 trial function is determined in 
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the same manner as N = 3, except that we now also in- 
corporate additional localization factors centered at sites 
{s 7 «/} which are situated near the six equivalent tertiary 
potential minima. This allows for a description of helium 
binding in these regions in the larger clusters. Three rep- 
resentatives of the sites s are marked on the 
potential contour plot of Fig. ^ The positions of all cor- 
responding sites are readily generated from application 
of the elements of the DQh(M) group. 

Cuts of the helium density along the y — plane are 
shown in Fig. These cuts are taken along a plane per- 
pendicular to the benzene molecular plane, bisecting two 
parallel carbon-carbon bonds. All structural quantities 
given here are obtained from DMC using the descendant 
weighting procedure described in Sec. lIIIBl and are thus 
sampled from the exact ground-state density distribution 
</>q. For N — 1,2, the helium density along the benzene 
Cg-axis has two maxima at z = ±3.7 A, which is in good 
agreement with the corresponding value of z = ±3.5 A 
inferred from spectroscopic measurements^ The posi- 
tions of these two density maxima remain unchanged as 
N increases to 14. As evident in Fig. [21 the local density 
near the benzene impurity is highly structured, which is 
reflective of the strong anisotropy in the helium-benzene 
interaction potential. We note that neglecting the rota- 
tional terms in the benzene kinetic energy, Eq. leads 
to a considerably more strongly peaked helium density 
distribution in the BF framed 



B. Excited states 

We have calculated helium vibrational excited-state 
energies for the HeAr-benzene system, using the POITSE 
methodology described in Sec. IIII CI For the N = 1 
dimer, each individual excitation operator A\ gives a 
well-defined excitation, by which we mean a spectral 
function n(E) consisting of a single, sharp peak. This 

indicates that the trial excited state given by Aj m ^T is 
a good approximation to a true eigenstate of the system, 
having non-negligible overlap with only a single <p n . The 
various spectral functions are superimposed and shown 
in the upper panel of Fig. [3] For N = 1, these helium vi- 
brational modes begin at ~ 9 K above the ground state. 
Below this onset, we would expect to see excitations in- 
volving rotational motion of the much heavier benzene 
impurity, which lie in the spectral range studied exper- 
imentally by Beck et alml The excitations shown in the 
upper panel of Fig. [3] can be grouped into pairs split by 
~ 0.6 — 2.3 K, corresponding to states generated by the 
application of pairs of projectors which are symmetric 
and antisymmetric with respect to reflection about the 
benzene molecular plane (Table ITTf> - Since the N = 1 
projectors used here all give well-defined excitations, we 
conclude that these states represent symmetric and an- 
tisymmetric tunneling pairs. We are not able to extract 
the lowest tunneling excitation given by the application 
of the i?io = z operator. This is likely due to its en- 



ergy being too close to the ground state Eq, and thus its 
DMC imaginary-time evolution is too slow relative to the 
DMC propagation times used for the excited-state calcu- 
lations here. These tunneling energetics are comparable 
to those obtained from basis set calculations of the 2,3- 
dimethylnaphthalene-He complex by Bach et al. 29 These 
authors found tunneling excitations associated with the 
motion of the single complexed helium atom from one 
side of the planar aromatic moiety to the other side, with 
splittings ranging from < 0.0014 K for strongly localized 
states, and up to 4.6 K for highly dclocalized states. 

The Bose-symmetrized versions of the excitation op- 
erators used for TV = 1 give a similar set of well-defined 
excitations for N = 2, which arc shown in the lower panel 
of Fig. |21 Thus, we conclude that these represent single- 
particle- like excitations, which is reasonable since the two 
helium atoms reside in the two equivalent global minima 
along the benzene C6-axis, and are well-separated by the 
benzene plane. In general, the N = 2 excitation energies 
E — Eq tend to be slightly lower than those for N = 1, 
with the onset of helium vibrational excitations begin- 
ning at ~ 7.5 K for N — 2, as compared to ~ 9 K for 
N = 1. Unlike the situation for N = 1, the addition of a 
second particle now allows for a well-defined excitation of 
symmetry at 7.63(4) K, obtained from the R\q = z 
operator. We note that neglecting the rotational terms 
in the benzene kinetic energy, Eq. J2J), alters the energy 
spectrum significantly^i 

With the two global minima occupied at N — 2, an ad- 
ditional third helium atom must be delocalized over the 
twelve secondary minima, due to the effect of helium- 
helium repulsions. Thus, the character of the excited 
states is expected to change dramatically as the num- 
ber of helium atoms increase from N — 2 to N = 3. 
This is evident in the POITSE calculations, where the 
projectors that give well-defined excitations for N = 1, 2 
(Fig. |3) now give multiple peaks for N = 3, shown in the 
upper panel of Fig. This indicates that the starting 
excited state ansatz generated by A, now nas appre- 
ciable overlap with more than one eigenstate. Since the 
POITSE methodology does not provide information on 
excited-state wave functions, interpretation of the N = 3 
excitations is now less straightforward. Additional in- 
sight can be gained by making modifications to the pro- 
jectors. The most noticeable new feature of the N = 3 
spectrum here is the appearance of lower-energy states 
at < 9 K. Since the two atoms situated near the global 
potential minima experience a very different local envi- 
ronment from the third atom that is delocalized over the 
twelve secondary minima, the question then arises as to 
whether one can ascribe any of the excitations to motion 
localized near the global potential minima alone. We 
explore this here by defining a set of weighted "local" 
excitation operators 

BLm = t u Rlmi "iy (52) 
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FIG. 2: Cuts of the helium density in the BF frame along the y — plane, for N — 1, 2, 3, 14. The two dominant peaks at 
z = ±3.7 A each correspond to a single helium atom situated on the molecule C6-axis, above and below the benzene plane. 



where the product in the denominator is taken over all 
localization factors that are not centered near the global 
minima. Application of these -B^-type projectors on ^>t 
gives the following local excited-state ansatz: 



(53) 



££(ri)£(r 2 ). 

3=1 



JV 



(54) 



By substituting Eq. l|4*5|) for £(r 3 -) m Eq. above, it 
can be seen that the action of Bj is to place an atom in 
an excited single-particle state that is spatially localized 
near the set of sites {s 7 <} only, i.e. near the two global 
potential minima. These operators are local in the sense 
that they act primarily on helium density near these two 
minima, while still maintaining spatial and Bose per- 
mutation symmetry. In contrast, by comparison with 
Eq. I]5ip. we see that the A' operators are global opera- 



tors, acting on all sites. 

The spectrum for TV = 3 generated by the set of W- 
type projectors is shown in the lower panel of Fig. 0] 
Compared to the spectrum derived from A^-type projec- 
tors (upper panel of Fig.^J, the spectral weight of lower- 
energy states are now reduced relative to the higher- 
energy excitations. A few low-energy features at < 9 K 
persist, in particular the R±q and R\\ excitations, albeit 
with reduced weights. Thus we conclude that the higher- 
energy excitations above ~ 9 K are associated with states 
that are spatially localized near the global minimum of 
the interaction potential, while the excitations below this 
range are associated with collective helium states of a 
more delocalized nature. For each set of states generated 
by a projector of given /, corresponding states of differ- 
ent m are clustered more closely together than in the 
N = 1, 2 spectra, with each cluster of different I spaced 
~3K apart. 

As N increases from 3 to 14, the ground-state DMC 
calculations reveal that the two global minima and twelve 
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FIG. 3: POITSE spectrum for TV = 1 (upper panel) and 
N — 2 (lower panel). The peaks are color-coded according 
to excitations derived from excitation operators of different I 
(Table EJ. 



secondary minima are completely occupied (see Sec. IV Al 
and Fig. |2J). An unambiguous determination of excited 
state energies at this size becomes more difficult due 
to the increase in Monte Carlo noise. However, cer- 
tain qualitative features remain apparent. The upper 
panel of Fig. [5] shows the POITSE excitation spectrum 
for N — 14, derived from the A'-type global projectors. 
These projectors generate trial excited states which over- 
lap with multiple eigenstates, and the resulting spectrum 
for a given projector has multiple peaks at both low and 
high energies. This problem, along with the accompany- 
ing increase in statistical noise in going to N = 14, gives 
a spectrum which is not fully converged in the sense that 
the peak widths are broader and the peak positions shift 
(by up to ~ 3 K) with additional sampling. Neverthe- 
less, we can again achieve qualitative insight by compar- 
ing with spectra derived from the localized Iv-type pro- 
jectors. The lower panel of Fig. [S] shows the spectrum 
derived from the B^-type projectors. There, the spec- 
tral weights of low-energy excitations (< 9 K) are again 
significantly reduced just as the N = 3 case, and each 
projector Bi m gives a k(E) consisting of a single domi- 
nant peak. While the higher-energy I = 4 excitations are 
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FIG. 4: POITSE spectrum for N = 3, deriving from it- 
type projectors (upper panel) and B -type projectors (lower 
panel). The peaks are color-coded according to excitations 
derived from excitation operators of different I (Table ITU. 



probably still not fully converged, the 1 = 1 — 3 peaks do 
not shift much (less than 0.5 K) with additional sampling, 
and appear to be converged. Thus the qualitative trends 
can be clearly established. Similar to the N — 3 case, 
the localized projectors give rise to a set of higher-energy 
excitations, except red-shifted by ~ 3 K in comparison 
to the N — 3 energies. For each I, states of different 
m are now clustered together as was also the case for 
N = 3 (Fig. 0} . These states represent the localized vi- 
brational motion of a single helium atom near the global 
potential minima, and appear to be associated with the 
localized helium density found in the path integral cal- 
culations of Ref. |2Sl The lower-energy excitations, on 
the other hand, are associated with collective vibrations 
delocalized around the periphery of the molecule surface. 



VI. SUMMARY DISCUSSION AND 
IMPLICATIONS FOR LARGER MOLECULES 

We have computed ground-state energies and struc- 
tures for small 4 He,/v-benzene clusters, where A?~ = 
1,2,3,14. For these sizes, the effect of the strong and 
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FIG. 5: POITSE spectrum for N = 14, deriving from it- 
type projectors (upper panel) and B -type projectors (lower 
panel). The peaks are color-coded according to excitations 
derived from excitation operators of different I (Table ITU. 



highly anisotropic helium-benzene interaction potential 
gives rise to a very structured helium density distribu- 
tion in the BF frame. In particular, a single helium atom 
is well-localized at each of the two equivalent global po- 
tential minima, above and below the benzene surface. 
We find a set of collective helium excitations with ener- 
gies of up to ~ 23 K above the ground state. Among 
these excitations, the higher-energy states (> 9 K) can 
be characterized as a localized excitation deriving from 
the helium density near the global potential minima, i.e. 
adsorbed on the molecular nanosurface. The existence of 
these localized modes is consistent with the localization 
of a single "dead" helium atom seen in Ref. |2f|. Helium 
excitations of lower energies were also obtained, which 
correspond to collective vibrations of helium atoms de- 
localized over equivalent sites of lower binding energy, 
situated near the edges of the benzene surface. Both 
localized and delocalized excitations can also be further 
classified by their symmetry with respect to the D§h (M) 
group. 

The energetic range of these helium excitations is simi- 
lar to that observed experimentally as vibronic structure 
in the mass-selective excitation spectra of planar aro- 



matic molecules in small helium clusters (N < 16)>i2iii 
It is also similar to the energetic range of the peaks ob- 
served in the phonon wings of vibronic spectra of planar 
aromatic molecules in large helium dropletsAL&S, Thus, 
these calculations provide support for the picture of he- 
lium atoms adsorbed and vibrating on the molecule sur- 
face, with the specific details of the helium motions being 
determined by the geometry of the surface. In particu- 
lar, the close relationship between the anisotropy of the 
molecule-helium interaction and the nature and energetic 
range of the excitations seen here for benzene implies that 
these molecule-induced vibrational modes will be very 
sensitive to the identity of the molecule, possibly even to 
the extent of providing a spectral "fingerprint" of com- 
plex polyaromatic species. Since these calculations have 
less than one solvation shell of helium surrounding the 
molecule, they are directly relevant to the recent exper- 
imental observations for small cluster sizesii2*ii More- 
over, the analysis of the excitations has been made in 
terms of single-particle type excitation operators acting 
on the many-body ground state, and is thus applicable 
to any number of solvating helium atoms. The fact that 
the localized excitation associated with the most strongly 
bound helium density at the global minimum persists in 
the largest cluster studied here (N = 14), together with 
the previous identification of localization at that site in 
larger helium clusters^ indicates that these localized 
excitations will also be present in much larger helium 
droplets. For the 4 HeAr-benzene system here we have ex- 
amined the limit of a single helium adsorbed atom on a 
single site given by the benzene nanosurface. But given 
the energetics that we find here, in a more general sense 
this class of surface-adsorbed vibrations are very likely to 
be responsible for the structure seen in the phonon wing 
sideband in experiments using helium droplets. 9 

A detailed analysis of the experimental phonon wing 
data clearly requires molecule-specific calculations with 
accurate potential surfaces. For example, the region of 
spatial confinement for an adsorbed helium atom near 
the global potential minimum is smaller in benzene than 
in the larger polyaromatic molecules studied experimen- 
tally with small numbers of helium atoms in Ref. In a 
general context, the present results are very encouraging 
in that they show that the diffusion Monte-Carlo based 
methodologies can be used to systematically study these 
excitations as a function of the number of helium atoms, 
provided that realistic molecule-helium interaction po- 
tentials are available. We emphasize here that a quanti- 
tative analysis, both for the small cluster vibronic exci- 
tations and for the phonon wing structure in electronic 
absorption spectra, will require detailed knowledge of the 
molecule-helium interaction potential in both ground and 
electronically excited states. 

Less obvious than the relation to phonon wing struc- 
ture is whether the excitations of the type studied here 
are responsible for the splittings on the order of ~ 1 K in 
the zero-phonon line, which were observed for some large 
organic impurity molecules in large helium droplets. This 
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zero-phonon line splitting has been experimentally stud- 
ied in detail for tetracene;^ and for indole derivatives^ 
In both experimental studies rotational fine structure has 
been ruled out as being responsible and several differ- 
ent interpretations have been advanced. For tetracene, a 
molecule possessing a plane of mirror symmetry, it has 
been suggested that either two inequivalent helium bind- 
ing sites exist, possibly due to an inhomogeneous solva- 
tion environment, or that some kind of two-particle tun- 
neling is manifested^. For indole derivatives, molecules 
possessing aromatic rings but no mirror symmetry and 
not usually completely planar, the theoretical evidence of 
localized helium atoms at aromatic rings in larger, super- 
fluid, helium clusters^ has been used to propose models 
of helium atoms similarly localized on either side of the 
aromatic portion of the molecules8> 

The 4 Hejv-benzene calculations reported here show 
low-lying excited states (< 9 K) for the larger clusters, 
N = 3, 14. Furthermore, the spectral weights of these 
low-lying states do decrease significantly relative to the 
higher-energy states when the local W-type operators 
are applied. This shows that the low-energy spectral fea- 
tures we observe the our calculations are due to vibra- 
tional motion of helium located near the periphery of the 
benzene surface, since the -B^-type operators specifically 
de-emphasize states that are not strongly localized near 
the global potential minima. In contrast to the higher- 
lying localized states, it is more difficult to say whether 
these low-energy states would remain low-energy with in- 
creasing N, since addition of more helium atoms would 
be expected to significantly change the local details of 
the helium wave function in these edge regions around 
the molecule. Thus, it is also more difficult to conclu- 
sively claim that these low-energy delocalized modes are 
responsible for the experimentally observed splittings in 
the zero-phonon lines. 

However, with aromatic molecules of lower overall sym- 
metry like the indole derivatives, it appears reasonable 
that the effective potential on opposite sides of the aro- 
matic ring may differ, due to the effect of non-symmetric, 
three-dimensional side chains, resulting in slightly differ- 
ent localized helium densities on either side of the aro- 
matic ring, as proposed in Ref. 0. This suggests that 
an alternative approach to interpret the splitting of zero- 
phonon lines in experiments would be in the context of 
impurity spectra in solids, where the impurity molecule 
can be trapped in structurally inequivalent trapping sites, 
giving rise to small differences in the spectral shifts of 



the electronic origin^ Quantitative analysis of this kind 
would also require determination of accurate helium sol- 
vation densities around the molecule. 

In extrapolating conclusions made from small cluster 
studies, whether theoretical or experimental, to explain 
experimental observations made in large droplets, the key 
point to consider is whether the localization induced by 
the helium-aromatic molecule interaction is sufficiently 
strong such that the localized helium excitations char- 
acterized here persist with increasing N. As the num- 
ber of helium atoms increases, collective compressional 
and surface modes of the droplet are expected to de- 
velop, and the coupling to these modes are manifest in 
the phonon-wing sideband in electronic spectra measured 
in large droplets^ The path integral results for benzene 
in 4 He39 provide evidence that the single localized he- 
lium atom located on cither side of the molecule near 
the global potential minimum does indeed retain its lo- 
calized identity in large clusters^ This provides a strong 
argument for the persistence of the single-particle- like lo- 
calized states whose excitations have been characterized 
here for smaller clusters (N < 14). An important task 
for future work is to understand how the discrete spectra 
for larger polyaromatics vary with the number of helium 
atoms when this is still less than a solvation shell^i as 
well as what happens to these excitations in much larger 
helium clusters. A theoretical prerequisite for study of 
all these excitations is now the development of accurate 
trial functions for these very anisotropic aromatic sys- 
tems at larger numbers of helium atoms (N > 14). This 
will then allow analysis of the transition from localized 
to collective helium vibrations, as well as identification 
of any persistent localized modes that may exist at an 
aromatic molecular nanosubstrate. 
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